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Cosmological Reionization Around the First Stars: Monte 
Carlo Radiative Transfer 
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ABSTRACT 

We study the evolution of ionization fronts around the first proto-galaxies by using 
high resolution numerical cosmological (A+CDM model) simulations and Monte Carlo 
radiative transfer methods. We present the numerical scheme in detail and show the 
results of test runs from which we conclude that the scheme is both fast and accurate. 
As an example of interesting cosmological application, we study the reionization pro- 
duced by a stellar source of total mass M = 2 x 1O 8 M turning on at z w 12, located 
at a node of the cosmic web. The study includes a Spectral Energy Distribution of a 
zero-metallicity stellar population, and two Initial Mass Functions (Salpeter/Larson). 
The expansion of the I-front is followed as it breaks out from the galaxy and it is 
channeled by the filaments into the voids, assuming, in a 2D representation, a charac- 
teristic butterfly shape. The ionization evolution is very well tracked by our scheme, as 
realized by the correct treatment of the channeling and shadowing effects due to over- 
densities. We confirm previous claims that both the shape of the IMF and the ionizing 
power metallicity dependence are important to correctly determine the reionization of 
the universe. 
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1 INTRODUCTION 

An increasing number of works has been recently dedicated 
to the study of the reionization of the universe (Gnedin & 
Ostriker 1997; Haiman & Loeb 1998; Valageas & Silk 1999; 
Ciardi et al. 2000, CFGJ; Miralda-Escude, Haehnelt & Rees 
2000; Chiu & Ostriker 2000; Bruscoli et al. 2000; Gnedin 
2000; Benson et al. 2000) which use both analytical and 
numerical approaches. An important refinement has been 
introduced by the proper treatment of a number of feed- 
back effects ranging from the mechanical energy injection 
to the H2 photodissociating radiation produced by massive 
stars (CFGJ). However, probably the major ingredient still 
lacking for a physically complete description of the reion- 
ization process is the correct treatment of the transfer of 
ionizing photons from their production site into the inter- 
galactic medium (IGM). Attempts based on different ap- 
proximated techniques have been proposed (Razoumov & 
Scott 1999; Abel, Norman & Madau 1999; Norman, Paschos 
& Abel 1998; Gnedin 2000; Umemura, Nakamoto & Susa 
1999), which sometimes are not readily implemented in cos- 
mological simulations. Therefore, it is crucial to develop ex- 
act and fast methods that can eventually yield a better treat- 



ment of the propagation of ionization fronts in the early 
universe. The first three papers are based on the ray-tracing 
method, whereas Gnedin (2000) uses the so called local op- 
tical depth approximation. Finally Umemura, Nakamoto & 
Susa (1999) implemented a time-independent ray-tracing 
method. Monte Carlo (MC) methods have been widely used 
in several physical/astrophysical areas to tackle radiative 
transfer problems (for a reference book see Cashwell & Ev- 
erett 1959) and they have been shown to result in fast and 
accurate schemes. Here we build up on previous experience 
of our group (Bianchi, Ferrara & Giovanardi 1996; Ferrara et 
al. 1996; Ferrara et al. 1999; Bianchi et al. 2000) in dealing 
with MC problems to present a case study of cosmological 
H11 regions produced by the first stellar sources. The study 
is intended as a test of the applicability of the adopted tech- 
niques in conjunction with cosmological simulations in terms 
of convergency, accuracy, and speed of the scheme. 

As an example of interesting cosmological application, 
we study the reionization produced by a stellar source of 
total mass M — 2 x 10 8 Mq turning on at z ~ 12, located 
at a node of the cosmic web. The study includes a Spectral 
Energy Distribution of a zero-metallicity stellar population, 
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and two Initial Mass Functions (Salpeter/Larson); the IGM 
spatial density distribution is deduced from high resolution 
cosmological simulations described below. 

In a forthcoming paper we will then improve the results 
of CFGJ by exploiting the strength of the MC method to 
release the approximations relative to the radiative transfer 
made in that paper. 



2 NUMERICAL SIMULATIONS 

We have studied the evolution of small scale cosmic structure 
in a standard A+CDM model with Q.a = 0.7 and Q. m — 0.3. 
The baryon contribution to fl m is Qt = 0.04, the adopted 
value of the Hubble constant is Ho = 70 km sec -1 Mpc -1 . 

The initial conditions are produced with the COS- 
MICS * code adopting a normalization ag = 1.14 for the 
power spectrum. The initial particle distribution (at red- 
shift z = 103.2) is made out of a 256 3 grid in a box of 
4 Mpc (comoving), thus giving a mass per particle of about 
1.5 x 1O 5 M0. We adopt vacuum boundary conditions; this 
requires the initial box to be reshaped. Assuming that tidal 
forces from large scale fluctuations are negligible for this 
problem, we cut a sphere of radius 2 Mpc from the initial 
cube and follow the evolution of the system in isolation. In 
practice, we study the evolution of a central sphere of ra- 
dius about 1.2 Mpc in full mass resolution with a coarse re- 
sampling of the most external particles. In the hi-resolution 
(central) region the number of particles is 2313847, whereas 
the rebinning left us with 10008 particles of variable mass 
(increasing with distance from the center) in the external 
1.2 - 2.0 Mpc shell. 

We evolve this initial matter distribution with GAD- 
GET (Springel, Yoshida & White 2000), a parallel N-body 
tree-code with an SPH scheme to describe hydrodynamical 
processes. SPH particles are placed on top of hi-resolution 
dark matter particles with the same velocity and a low ini- 
tial temperature (about 150 K, the temperature of inter- 
galactic medium as determined by adiabatic expansion of 
the universe after decoupling); the mass of each fluid parti- 
cle is rescaled according to the baryon and matter fraction 
quoted above. The (spline) softening length is 1 kpc for the 
gas particles, and 3 kpc for hi-resolution dark matter par- 
ticles, respectively. For the problem at study here, the sim- 
ulation is stopped at z ~ 5, thus avoiding the very strong 
nonlinear phase in which our assumption of zero tidal effects 
from large scale structure would break down. 

The gas-dynamics evolution is purely adiabatic in order 
to reasonably limit the integration time. Thus, we are not 
able to follow the details of cooling processes that ultimately 
lead to the formation of very dense and cold star forming re- 
gions. However, the density field of the mildy overdense IGM 
should be reasonably well described by this assumption. Our 
mass resolution does not also allow a full resolution of the 
Jeans mass, and therefore the presence of sub-grid structure 
that we cannot track is expected. This is a problem com- 
mon to all presently available simulations, implying that the 
treatment of IGM clumping is only approximate, as recently 
pointed out by Haiman, Abel & Madau (2000). 

* Bertschinger, http://arcturus.mit.edu/cosmics/ 
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Figure 1. Evolution of the SED for a Simple Stellar Population of 
metal free stars for two different IMFs: Salpeter (solid lines) and 
Larson with M c = 5Mq (dotted). The evolutionary times are, 
from the top to the bottom, 5, 7, 10, 15, 30, 50, 70, 100, 150 Myr. 



Simulations were performed on the Cray-T3E at 
Rechenzentrum-Garching, the Joint Computing Center of 
the Max-Planck-Gesellschaft and Max-Planck-Institut fuer 
Plasmaphysik. 

Finally, from the SPH particle temperature and density 
distribution we reconstruct a central cubic region of 128 3 
cells, the grid values of any field being evaluated from par- 
ticle values with an SPH-like interpolation. For the reion- 
ization study presented below, we have concentrated, as an 
example, on a source located at redshift z « 12 immersed in 
the density field calculated from the simulations described 
above. The location and mass of the dark matter halos in 
the simulation have been determined by means of a friend- 
of-friend algorithm. We have then selected a halo suitably 
located in the node of a cosmic filamentary structure and 
associated a luminous source to it. The total halo mass is 
M = 2 x 10 s Mq\ the corresponding stellar mass has been 
determined following the same prescriptions of CFGJ, i.e. 

AL = (lV«m)M/i,/*, (1) 

where fb = 0.08 is the fraction of virialized baryons able to 
cool and become available to form stars, and /* = 0.15 is the 
star formation efficiency. Although with our cell size of 15 
kpc we are marginally able to resolve the source (comoving) 
virial radius (about 20 kpc) , we do not attempt to treat the 
the escape of ionizing photons (Dove, Shull & Ferrara 2000; 
Wood & Loeb 2000) from the galaxy interstellar medium. 
Instead we take a fixed value for this quantity equal to / esc = 
0.2, the upper limit obtained by the above studies, by which 
we multiply the source ionizing photon flux. 
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3 POPIII STAR SPECTRUM AND IMF 

The SEDs of metal free star systems used in this paper were 
calculated implementing the synthetic evolutionary code de- 
veloped for a Simple Stellar Population (Brocato el al. 1999, 
Brocato el al. 2000). For a detailed description of the com- 
putational procedure we refer to those papers; here we just 
outline that the results of this code have been largely tested 
on young LMC clusters and old galactic globular clusters 
showing a very good agreement. The present synthesis mod- 
els rely on the set of homogeneous evolutionary computa- 
tions for stars with a cosmological amount of heavy elements 
and helium (Z=10~ 10 and Y=0.23) presented by Cassisi & 
Castellani (1993) and Cassisi, Castellani & Tornambe (1996) 
and properly extended to higher masses in order to cover a 
wide range of ages. All the major evolutionary phases are 
taken into account, both the H and the He burning phase 
for stars with original masses in the range 0.7—40 Mq, for 
which the He burning phase is followed up to the Carbon 
ignition or, alternatively, to the onset of the thermal pulses. 
Models atmospheres are by Kurucz (1993) for metallicity 
Z=2xl0 -7 . The available temperatures range from 3700 K 
to 50000 K, so the spectra at higher temperatures and grav- 
ity are simply extrapolated in order to keep the homogeneity 
of the grid. The wavelength range is from 91 A to 160 
and the spectral resolution is typically 10-20 A in the UV 
to optical bands. We consider an instantaneous burst of star 
formation assuming that masses range from 1 Mq, as re- 
cently suggested for Pop III stars by Nakamura & Umemura 
(1999a, 1999b), to 40 M and two cases for the IMF: i) a 
standard Salpeter law, and ii) a Salpeter function at the 
upper mass end which falls off exponentially below a char- 
acteristic stellar mass M c (Larson 1998). The last case is 
taken into account since both theory of thermo-dynamical 
condition of primordial gas and the absence of observations 
of solar mass, metal free stars seem to suggest at birth a 
paucity of low mass stars. Moreover, the mass interval used 
in this paper properly describes population models of age 
down to 5 Myr, the lifetime of a 40 Mq star; therefore, the 
burst models of this age are actually populated by 1-40 Mq 
stars. The resulting grid of SEDs were calculated at intervals 
of 1, 10 and 50 Myr between 1 and 20 Myr, 20 and 100 Myr, 
and 100 and 150 Myr, respectively. Other calculations of 
zero-metallicity star SEDs are present in the literature (see 
e.g. Cojazzi el al. 2000; Tumlinson & Shull 2000), but they 
do not show the detailed SED time evolution. 

In Fig. 1 we show the adopted SEDs for the Salpeter 
IMF and the Larson IMF with M c = 5M . The two sets of 
curves do not differ sensibly in their spectral shape evolution 
but the Larson IMF has about 4 times larger specific photon 
flux at the Lyman continuum. In Fig. 2 we show the time 
evolution of the total number of ionizing photons (in units 
of 10 60 ) emitted by a solar mass of metal free stars with 
the two IMFs. After about 10 7 yr both curves flatten as a 
consequence of the decreasing number of surviving massive 
stars. It is important to note that, during the first 5 Myr, a 
zero-metallicity stellar cluster with a Salpeter IMF has an 
ionizing photon rate about 25 times higher than the anal- 
ogous one for a 1% solar metallicity used by CFG J. As we 
will see later, this has important effects on the reionization 
process. 
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Figure 2. Time evolution of the total number of ionizing photons 
(in units of 10 60 ) emitted by a solar mass of metal free stars 
with a Salpeter (solid line) IMF and a Larson (dotted) IMF with 
M c = 5M Q . 



4 MONTE CARLO RADIATIVE TRANSFER 

The application of MC schemes to radiative transfer prob- 
lems requires that the radiation field is discretized in a rep- 
resentative number of photon packets, A/j,. The processes 
involved (e.g. packet emission and absorption) are then 
treated statistically by randomly sampling the appropriate 
distribution function. Let us consider a source of bolometric 
luminosity L and lifetime t s (L can be a function of time 
as our method easily allows to treat short-lived and vari- 
able sources). Packets have the same energy, £ p , but contain 
a different number of monochromatic photons, 7V 7 , of fre- 
quency v so that for the j-th packet it is E p = N~{(j)hv(j); 
their rate is AT P — L/£ p and the simulation time when the 
j-th photon packet is emitted is then I = jdt (j = 1, .., Af p ), 
where dt — AC" . The number of photons in each packet 



emitted by the source during t s is: 



E s 



Lt a 



N p hv(j) N p hv(j) ' 



(2) 



The packet frequency, v, is obtained by sampling the source 
SED, S(y): if R v is a random number, v is given implicitly 
by: 



dv'Siy') 



dv'Siy') 



(3) 



where vh is the threshold for hydrogen ionization. 

We adopt spherical coordinates (r, 9, </>) with origin at 
the source location (x s ,y s ,z s ) in the box, such that the 
cartesian coordinates are: 



x = r sin# cosfi 0^0^ 2-7T, 
y = rsin#sin<^ ^ 9 ^ n, 
z — r cos#. 



(4) 
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Assuming that the emission is isotropic and spherically 
symmetric, photon packets will propagate along the direc- 
tion: 

s = [r, arccos(l - 2R e )6, 2nR t/> <f>], (5) 

where Re and R^ are random numbers. Note that the above 
choice automatically ensures a constant surface density of 
photon packets, thus preventing flux anisotropics at the 
poles due to angle sampling; this can be seen as follows. As 
we are sampling Re from a uniform distribution, the fraction 
of photon packets emitted in the angular range 8\ < 8 < 82 
is equal to 

dM p (6 1 ,8 2 ) = {Re 2 -Re 1 ), (6) 

where Re i = (1/2)(1 — cos 9i) from eq. 5. As the area fraction 
of the sphere belt delimited by 81 and 82 is 

dA(0i,0 2 ) = ~(cos0i-cos0 2 ) (7) 

one sees that the surface density dN v /dA of photon packets 
is constant and independent of angles. 

Once emitted, a packet of frequency v will travel for a 
finite path length of optical depth: 

i i 

7-iO) = a{v) ^(iVHi) fc = e(v) j^(T»m)*/Ax, (8) 

k=l h=l 

before being absorbed in the i-th grid cell. Here a(y) is 
the photoionization cross section, (Nm)k — (nm)k.fAx is 
the IGM neutral hydrogen column density through the k- 
th cell of size Ax; the cell with k = contains the source 
location and it is not included in the opacity count, as 
its contribution is already included in f esc . For simplicity, 
the results presented here only include H opacity; inclu- 
sion of He, molecules (as H2 and HD) and heavy elements, 
if present, is straightforward, although computationally ex- 
pensive for molecules as one has to properly compute both 
the line widths and their Doppler shifts. The factor / ac- 
counts for the fact that the path through a cell is in the 
range < £ ^ \^3Ax depending on its inclination. To limit 
the computational cost, we do not attempt to calculate the 
actual path and we treat this effect statistically as follows. 
Given a cubic box of unit size Ax — 1, via an indepen- 
dent Monte Carlo procedure, we randomly pick up on a face 
the coordinates of an entering ray and the direction of its 
propagation. We then derive the differential probability dis- 
tribution function, P(£), of the path lengths (Fig. 3) and 
define / = 0.56 as the median value of such distribution. 

The value of i in eq. 8 is defined as the minimum one 
for which the following condition is satisfied: 

n(y) > T = -ln(l - Rr), (9) 

with R T a random number. The previous relation is obtained 
by sampling the photon absorption probability distribution. 

In each cell k — 1, along the path from the source 
to the absorption site of the j-th packet, we update the 
value of the hydrogen ionization fraction; x = n H + /tih, %+ 
and riH are the ionized and total hydrogen number density, 
respectively. This is done by advancing the time-dependent 
ionization equation in discretized form. At time t = t n+1 it 
is: 

xt +1 = x n k + AtMxt) - 7r (a£)] + N ^ 5 ki ; (10) 
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Figure 3. Differential probability distribution function of 
lengths, £ for randomly oriented paths entering from a face of 
a cubic box of unit size. The integral of the curve is normalized 
to unity. 

7 C and y r are the collisional ionization and recombination 
rates, respectively (Cen 1992); they are evaluated at the 
temperature of the gas in the cell. The initial temperature 
is increased to T = 10 4 K after absorption of a photon 
packet and set equal to T = 2.734(1 + z), i.e. the CMB 
temperature, after the cell recombines (see below). This is 
a reasonable approximation, although T could be smaller 
at higher redshift (Ciardi, Ferrara & Abel 2000). The last 
term, describing photoionization due to photon packet ab- 
sorption, is present only for the cell in which absorption 
takes place (k = i, Ski = 1); in the other cells along the path 
(k < i, Ski = 0) we simply follow the gas recombination. In 
case the derived value of x^ +1 > 1, we set x^ +1 — 1, and 
we use the extra photons to ionize the cell Xk+i- In practice, 
we write Nj(j) = N$fi(j) + N~(j), where the number of 
photons required to obtain x^ +1 = 1 is: 

N( 1) (j) = {(l-x%)-Ath c (xt)- 7r (xl)]}(n ii ) k (Ax) 3 .(ll) 

The most delicate point of the ionization fraction evaluation 
is to find an appropriate expression for At. In principle, one 
could update the ionization state on the entire grid after 
each photon packet emission. However, this is computation- 
ally too expensive. For this reason we choose to update the 
ionization state of cells along the packet path and we deduce 
At as follows. 

The cell in which the source is located is enclosed by a 
cubic surface (shell). If we define the order of this shell to 
be g = 1, the number of its member cells is given by: 

Af s ( 5 ) = (2 5 + l) 3 -(2 5 -l) 3 = 2(12 3 2 + l) g>l, (12) 

or N s (l) = 26. The second order shell (g = 2) is instead 
made of 98 cells and so on. As the photon packet directions 
are isotropically distributed, statistically a cell in a shell of 
order g will be along a path every [N a (g)]~ 1 packets. There- 
fore the typical time scale between two subsequent ionization 
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fraction updates is At — N s (g)dt; this is the value we use in 
eq. 10 above. The order g is calculated as: 



g — Max[k x 



(13) 



where the fc;'s indicate the discrete cartesian coordinates of 
the cell in the reference frame eq. 4. The above approxi- 
mation sets a lower limit on M p , as in order to be valid it 
is necessary that the update time At is shorter than the 
average recombination time, (t r ), in the box. Hence to ad- 
vance the simulation up to a time t„ the number of packets 
required is: 

AA p >24iV 2 if y , (14) 

where iV 3 is the total numer of cells in the box. This detailed 
treatment of recombination is important to accurately eval- 
uate the optical depth through eq. 8. 

We then introduce an ionization vector (size N 3 ) whose 
elements are initially set to zero. The vector element cor- 
responding to a cell for which the ionization fraction is in- 
creased above x t h = 0.99 as a result of an absorption event 
is updated to contain the value: 



Jr 



Int 



t r {x k , n k ) 
dt 



(15) 



where t r (xk,Tlk) is the local recombination time and t is 
the current simulation time. If further absorption does not 
take place before the time j r dt, we allow for recombination 
emission from the given cell; if instead another absorption 
event takes place, we re-update the vector element to a new 
value of j r . As already mentioned, after recombination, the 
cell temperature is set at the CMB temperature at that red- 
shift. 

Following a cell recombination a packet is emitted along 
a random direction from the same location. The probabil- 
ity that the emitted packet has hu(j) > 13.6 eV is given 
by the ratio between the number of recombinations to the 
ground level and the total one, i.e. 50% for T = 10 4 K. The 
frequency of the re-emitted packet is determined through 
eq. 3, where the source SED is substituted by the volume 
emissivity of the diffuse radiation, given by the Milne rela- 
tion: 



S d (u) = 



2hv i 



3/2 



g i+1 \ v 27rm p fcT ( 

, \ -h(u-u H )l(kT B ) 

<j(v)e 1 



n H +n e 



(16) 



where gi and gi+i are the hydrogen ground level statistical 
weights; T e and n e are the electron temperature and num- 
ber density. The propagation and absorption of re-emitted 
packets is then followed in the same manner as for primary 
ones. The only difference is that the opacity contribution of 
the re-emitting cell is now included. 

The method is easily generalized to include an arbitrary 
number of point sources once the packet emission rates, and 
hence their emission sequence for the various sources are 
given. 

4.1 Testing the Method 

The method described in the previous Section has been 
tested against the analytical solution for the temporal evo- 
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Figure 4. Residuals of the numerical Hji region equivalent ra- 
dius, R n , for the convergency test: crosses refer to D 5 , circles to 
D 50 (see text for definition of D). The source has ionizing photon 
rate J\f~ = 2.5 X 10 50 s -1 . 



lution of the radius of an Hn region produced by a source 
with constant ionizing rate and expanding in an homoge- 
neous medium (Donhaue & Shull 1987; Shapiro & Giroux 
1987). The Hn region equivalent radius, R n , is numerically 
derived by first assuming that a cell k is in the ionized state 
when Xk > Xth = 0.99. If Ni on is the number of ionized 
cells in the box, the volume occupied by the ionized re- 
gion is V = Ax 3 Ni 0n and the equivalent radius is then 
i?n = (3y/47r) 1 ^ 3 . In all the tests presented here the gas 
density is set equal to average IGM one in an Einstein- 
de Sitter universe with Q,^h 2 = 0.019 at z = 11.6, i.e. 
nu = 4.5 x 10~ 4 cm -3 ; the adopted spectrum is monochro- 
matic, with a photon energy equal to 13.6 eV. 

First, we check for convergency of our scheme. We run 
three simulations with different values of M p = (5, 50, 500) x 
10 5 . The source ionizing photon rate is M 1 = 2.5 x 10 50 s _1 . 
For each run we calculate the residuals with respect to the 
highest resolution run M p = 500 x 10 ; they are defined as: 



D y 



R v - R 5 

p500 
1 



(17) 



where the superscript y — 5, 50 refers to the low- and 
medium-resolution runs, respectively. The residuals are 
shown in Fig. 4 at different evolutionary times together 
with the associated mean quadratic errors; the error on 
R n is assumed to be equal to Ax. Satisfactory conver- 
gency is achieved by the medium-res run; subsequent tests 
have hence been done using Af p = 50 x 10 5 . For reference, 
this run takes only ~ 300 s (CPU time) on a SUN UL- 
TRA10/333MHz workstation. 

Next, we assess the accuracy of the solution. In Fig. 5a 
we show the time evolution of R n , as produced by sources 
of different ionizing photon rates A/" 7 = (2.5, 12.5, 25, 125) x 
10 50 s _1 , compared to the corresponding analytical radius, 
R a . Fig. 5b shows the residuals (defined analogously to 
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Figure 5. (a) Time evolution of the equivalent Hn radius 
produced by sources of different ionizing photon rate J\fy = 
2.5(circles), 12.5(squares), 25(triangles), 125(crosses) X 10 50 s" 1 , 
compared to the analytical solutions (solid lines), (b) Residuals 
of the numerical radii calculated with respect to the analytical 
ones. The notation is the same as in Fig. (a). 



eq. 17 with substituted by R a ) with respect to the an- 
alytical solution. In the worst case (very faint source, early 
evolutionary times) the residual is about 20%; This rela- 
tively large discrepancy occurs because the size of the ion- 
ized region is only a few cells. At later times and/or for 
more luminous sources the agreement is remarkably good, 
i.e. within a few percent. 

Finally, we check the effects of varying the spatial res- 
olution. To this aim we have considered runs with dif- 
ferent box sizes, keeping a fixed linear number of cells 
N — 128. In Fig. 6 we show the residuals for runs with 
a box size (1,0.2,0.1) x NAx with respect to the high- 
est spatial resolution run for a source ionizing photon rate 
ACy = 2.5 x 10 50 s _1 . 
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Figure 6. Residuals for runs with a comoving box size (1, 0.2) X 
NAx (circles and crosses, respectively) with respect to one with 
N Ax/10. The source ionizing photon rate is Afy = 2.5 X 10 50 s" 1 . 



5 RESULTS AND DISCUSSION 

The combined cosmological and radiative transfer simula- 
tions described above allow us to determine the evolution of 
the ionized region around the selected zero-metallicity stel- 
lar object. Illustrative slices extracted from our simulation 
box through the source location are shown in Fig. 7 at dif- 
ferent times (10, 20, 40, 60, 100 Myr) after the source has 
been turned on and a for Larson IMF, together with the ini- 
tial density field. The above time interval corresponds to less 
than two redshift units at z ~ 12, and it is therefore not un- 
reasonable for the purposes of this paper to neglect the IGM 
density evolution As seen from Fig. 7, the ionization front 
(I-front) breaks out from the galaxy very rapidly, leaving 
behind a very clumpy ionization structure in the immediate 
surroundings where the IGM is overdense and, consequently, 
recombination times are shorter. As an aside, it is interesting 
to note that, as the circular velocity of the parent galaxy is 
a few km s _1 , the ionized gas (whose sound speed is of order 
of 10 km s~ ) will very likely be able to leave the galaxy, 
thus quenching further star formation. If this is a widespread 
phenomenon, it might have strong implications for the evo- 
lution of dwarf galaxies as already outlined by some authors 
(Barkana & Loeb 1999; Ferrara & Tolstoy 2000). Once the 
I-front expands in the IGM, it is channeled - for similar rea- 
sons - by the large filaments into the underdense volumes 
(voids) and its structure becomes very complex and jagged, 
assuming, in a 2D representation, a characteristic butterfly 
shape. The typical overdensity of the clumps encountered 
by the expanding I-front is « 30. Qualitatively, the typi- 
cal final extent of the ionized region is of about 1 comoving 
Mpc; this size is reached already after 60 Myr. After that 
time, the rapid decrease of the source ionization power (Fig. 
2) slows down the expansion. The final stage of the evolu- 
tion is constituted by a relic Hn region which slowly starts 
to recombine on a time scale which in the voids can be as 
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Figure 7. Slices through the simulation box containing a metal 
free stellar source of total mass M = 2 X 10 s Mq at z ss 12 in 
a A+CDM cosmological density field and a Larson IMF. The six 
panels show the neutral hydrogen number density, superposed to 
which (black regions) is the evolution of the ionized (x > = 
0.99) hydrogen distribution at times 0, 10, 20, 40, 60, and 100 
Myr after source turn on. The linear (proper) size of each slice is 
160 kpc. 



long as 0.8 Gyr. The MC technique shows here all his power 
in following the details of the I-front evolution. For exam- 
ple, it tracks remarkably well the channeling induced by the 
large scale structure: the ionization front cannot propagate 
inside the densest filaments due to their large optical depth 
and short recombination time. In addition, also the effects 
of shadowing produced by isolated clumps are clearly recog- 
nized from the fingers protruding from the Hn region to the 
left of the source in the upper panel of Fig. 8. The ionization 
cone visible below the source, caused by a large overdensity 
located close to the source in that direction, is also a result 
of shadowing. Thus it seems that our scheme is highly suit- 
able at least for this type of cosmological radiative transfer 
calculations. 

We now turn to the differences induced by our two as- 
sumptions concerning the IMF. Fig. 8 shows a comparison 
between the final stages (100 Myr after the source turn on) 
of the I-front evolution for a Larson (upper panel) and a 




Figure 8. Comparison between the final stages [t = 100 Myr) of 
the I-front evolution for a Larson (upper panel) and a Salpeter 
(bottom) IMF for the same source properties as in Fig. 7. 



Salpeter (bottom) IMF. The source stellar mass and other 
properties of the simulations are the same as those discussed 
above. For a Salpeter IMF, the volume of the ionized region 
is smaller by a factor 8, although the shape is very similar to 
the one for a Larson IMF case at an earlier stage, roughly 
corresponding to 10 Myr (see Fig. 7). This was expected 
from the two adopted SEDs. In fact, the total number of 
ionizing photons (see Fig. 2) per stellar mass formed inte- 
grated over the entire source lifetime and spectral extent is 

5 x 10 B1 (10 61 ) for the Larson (Salpeter) IMF. Thus, the 
IMF might play an important role for the reionization of 
the universe; in addition, zero-metallicity stars have larger 
ionizing power as already stressed previously and recently 
addressed by other authors (Cojazzi et al. 2000; Tumlinson 

6 Shull 2000). 

To assess the differences between the evolution of an 
I-front propagating in an inhomogeneous medium, as in 
the simulations discussed here, and in an homogeneous gas 
at the mean IGM density, we have calculated the ratio 
V = V/Vh between the ionized volumes in our simulation 
and in the homogeneous case. This is found to be roughly 
independent of time and equal to V = 0.45 for the Larson 
IMF. In the Salpeter IMF case V varies slightly with time 
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from 0.125 to 0.2. Hence, this result confirms that the ion- 
ized volume tends to be smaller in the inhomogeneous case, 
but suggests that the effect is not dramatic. 



ACKNOWLEDGMENTS 

We would like to thank S. Cassisi for providing us with 
the stellar tracks; N. Yoshida and V. Springel for advises 
in various numerical problems; the referee T. Abel for use- 
ful comments. This research was supported in part by the 
National Science Foundation under Grant No. PHY94-07194 
(SM) and by the Italian Ministery of University, Scientific 
Research and Technology (MURST) (GR). 



Shapiro, P. R. & Giroux, M. L. 1987, ApJ, 321, L107 

Springel, V., Yoshida, N. & White, S. D. M. 2000, preprint (astro- 
ph/0003162) 

Tumlinson, J. & Shull, J. M. 2000, ApJ, 528, 65 

Umemura, M., Nakamoto, T. & Susa, H. 1999, in Numerical As- 
trophysics, eds. Miyama et al. (Kluwcr: Dordrecht), p. 43 

Valageas, P. & Silk, J. 1999, A&A, 347, 1 

Wood, K. & Locb, A. 2000, ApJ, in press 



REFERENCES 

Abel, T., Norman, M. L. & Madau, P. 1999, ApJ, 523, 66 
Barkana, R. & Loeb, A. 1999, ApJ, 523, 54 

Benson, A. J., Nusscr, A., Sugiyama, S. & Lacey, C. G. 2000, 

preprint (astro-ph/0002457) 
Bianchi, S., Ferrara, A. & Giovanardi, C. 1996, ApJ, 465, 127 
Bianchi, S., Ferrara, A., Davies, J. & Alton, P. 2000, MNRAS, 

311, 601 

Brocato, E., Castcllani, V., Raimondo, G. & Romanicllo, M. 1999, 
A&AS, 136, 65 

Brocato, E., Castellani, V., Poli F. M. & Raimondo, G. 2000, 

A&A, submitted 
Bruscoli, M., Ferrara, A., Fabbri, R. & Ciardi, B. 2000, MNRAS, 

318, 1068 

Cashwell, E.D. & Everett, C.J. 1959, A Practical Manual on the 
Monte Carlo Method for Random Walk Problems (New York: 
Pcrgamon) 

Cassisi, S. & Castellani, V. 1993, ApJSS 88, 509 

Cassisi, S., Castellani, V. & Tornambe A. 1996, A&A, 317, 108 

Ccn, R. 1992, ApJS, 78, 341 

Chiu, W. A. & Ostriker, J. P. 2000, preprint (astro-ph/9907220) 
Ciardi, B., Ferrara, A. & Abel, T. 2000, ApJ, 533, 594 
Ciardi, B., Ferrara, A., Governato, F. & Jenkins, A. 2000, MN- 
RAS, 314, 611, CFGJ 
Cojazzi, P., Bressan, P., Lucchin, F., Pantano, O. & Chavez, M. 

2000, MNRAS, 315, 51 
Donahue, M. & Shull, M. J. 1987, ApJ, 232, L13 
Dove, J. B. , Shull, J. M. & Ferrara, A. 2000, ApJ, 531, 846 
Ferrara, A., Bianchi, S., Dettmar, R.-J. & Giovanardi, C. 1996, 
ApJL, 467, 69 

Ferrara, A., Bianchi, S. Cimatti, A. & Giovanardi C. 1999, ApJSS, 
123, 437 

Ferrara, A. & Tolstoy, E. 2000, MNRAS, 313, 291 
Gnedin, N. Y. 2000, preprint (astro-ph/0002151) 
Gnedin, N. Y. & Ostriker, J. P. 1997, ApJ, 486, 581 
Haiman, Z. & Loeb, A. 1998, ApJ, 503, 505 

Haiman, Z., Abel, T. & Madau, P. 2000, preprint (astro- 

ph/0009125) 
Kurucz, R. L. 1993, CD Rom n. 13 
Larson, R. B. 1998, MNRAS 301, 569 

Miralda-Escudc, J., Haehnelt, M. & Rees, M. R. 2000, ApJ, 530, 
1 

Nakamura, F. & Umemura, M., 1999a, A&A, 343, 41 
Nakamura, F. & Umemura, M., in Star Formation 1999b, Proc. 

of Star Formation 1999, held in Nagoya, Japan, Ed. T. 

Nakamoto, Nobeyama Radio Observatory, p. 28-29 
Norman, M. L., Paschos, P. & Abel, T. 1998, in H 2 in the Early 

Universe, Florence, Proc., eds. E. Corbelli, D. Galli & F. Palla, 

p. 455. 

Razoumov, A. & Scott, D. 1999, MNRAS, 309, 287 



© 2000 RAS, MNRAS 000, 1-8 



